|
Spring Mass Damper Simulation Tutorial
Springs usually occur physically as a coil of metal, and their idealizations have pretty simple behavior: compressing the spring will result in the spring pushing back, and stretching the spring will have it trying to pull back towards the start position, so any displacement along the axis of the spring will be countered by an opposite force that will tend to move the spring back to it's original position. Any basic mechanics/physics book will tell you the fundamental spring equation: F = -kxWhere k is the spring constant (how loose or springy the spring is), x is the difference between the springs current length and its rest length, and F is the force on both endpoints of the spring. Usually one endpoint is fixed, the other is the one that bounces around- which is usually what happens: an initial impulse displaces the spring, the unfixed end of the spring acquires some velocity moving back, but it passes through the zero-displacement point, is pulled back in the other direction, and may bounce forever in the absence of any dampening forces. Physical springs have more complex behavior(like the transverse vibration and accompanying doing! sound when they're bent away from their axis) and could be described by more complex models but we'll start from the simplest model. The first step is to implement the spring equation in 3D. F and x in the equation are vectors, and k is a scalar. A well designed vector class will take care of most of the 3D details for us, so we can find x with something like:
vector spring = end[1]->location - end[0]->location;
Where we're inside our spring class, the spring has an array of two pointers to ends (masses, actually, which may be connected to other springs) with location vectors. We then take the length of this displacement vector (l = sqrt(x^2 + y^2 + z^2) and subtract the original length:
float length = spring.Length();
Since displacement is a scalar here and we'll need to orient the force along the spring axis we can get the normal:
vector springN = spring/length;
(Again, the vector class has been overloaded to accept vector/scalar operations) Now comes the fundamental spring equation:
restoreForce = springVecN*((displacement*K);
This force is imparted on both ends of the springs, but in opposite directions. With a little thought or experiment (the spring will simply expand and never bounce back if you get it wrong) you can get the right directions:
end[0]->force = end[0]->force + restoreForce;
We'll update the end masses every cycle with:
velocity += force/mass;
(Also remember to zero out the force after the update, lest the springs keep accumulating force without bound) That's it for the simplest model, and the easiest way to test is to hook up two masses to one spring, and fix one of the masses (don't update its position). Impart some impulse on free end and then watch what happens...
The simple model seems bounce correctly, but then something odd happens: when the spring completes one cycle, it extends further out than the initial impulse pushed it, and then has that much more energy that it bounces back even further out in the next cycle. It goes out of control, oscillating uncontrollably and violating the law of energy conservation(which doesn't exist inside our system). What went wrong? There's a couple of things that make this more complicated than I originally proposed: We update our system a finite number of times during a spring oscillation, while the physical world 'updates' continously. Each updates each impart forces on the mass that are then tell the masses how to change velocity for one time step, but really any arbitrarily small movement should result in a recomputation of spring forces, and new forces will tell our masses how to move differently. We can improve upon this by increasing the number of updates, but our computing resources and desire for other processing intensive happenings will limit how much we can do this. We could also solve the equations (with derivatives and integrals and so forth) that would tell us what forces are going to exist in the time inbetween updates, but this gets more difficult (if not impossible) for every additional attached spring. Probably our floating point coordinates that make up the vectors are only 32 bits. Operations that result in trailing digits are just rounded or cut off, and the more operations the more errors that will be accumulated. This is probably a much smaller problem than the first, but should be considered. There's a method call Range-Kutta that could does a sort of psuedo differentiation that might benefit the system, but for now we'll move on to dampers.
A damper is kind of the opposite of a spring, except it operates on relative velocity rather than displacement. Spring endpoints moving away from each other will have forces imparted from the damper that will act against that motion (only on the spring axis, however), as well as endpoint moving towards each other. This will tend to return the spring to a static position. Also endpoints moving in unison will not be affected (the damper won't act as drag), and one endpoint unmoving and the other moving will average out to both moving slower than the one endpoint. Conceptually, we first need to find the velocities along the spring axis: Click to view imageAnd then what behavior do we want our damper to exhibit? Click to view imageI found this equation in the O'Reilly book "Physics for Game Developers", which didn't seem to mention the axis alignment issues but basically fits requirements of the table above: Fd = (V1-V2)*kdIn the actual code we can find the dot products after finding the change in velocity, which will be equivalent to taking the dot's then finding delta V.
deltaV = end[1]->velocity - end[0]->velocity;
We have a damper constand Kd to work with:
damperForce *= Kd;
d is a scalar, so to we multiply by the spring normal again to get direction:
damperForce = springVecN*d;
A composite equation looks like:
damperForce = springN*((diff*K) + Dot(&springN,&deltaV)*Kd);
Now will it work?
The most time-consuming thing to do now is to create objects that will test the system. The single spring system should also be tested with a constant force (rather than a single endpoint), where the spring will bounce and then find a stable point where 'gravity' is matched by the spring restoring force. Also, if gravity is perpendicular to the initial spring orientation, the free mass should swing back and forth like a pendulum. Hook up a 3d shape and see if perturbations make the shape collapse or if the setup is stable and the shape restores itself. A critical requirement for all of this is that behavior is exactly the same for when the system has some average velocity as when the average velocity is zero- which suggests testing the system by do just that (perhaps with a moving camera that will keep the system on screen, and a stationary background that will give the impression of motion). If you perturb a pyramidal (for example) shape just right it will begin to rotate, and problems arise: we aren't doing any rotational dampening, so roundoff errors may start the spring rotating increasingly faster (and the springs will stretch out as if we had implemented centrifugal 'force', which isn't a force at all in the cartesian system and which is why it doesn't need to be explicitly implemented in order to arise) and the pyramid blows up. A simplistic 'drag' way to fix this is to just dampen all velocity:
velocity *= .99;
It would be nicer if our system could somehow recognize instability situations and then dampen accordingly (as the above dampening will dampen linear motion as it does rotational, which is bad if you're trying to model bodies in space that shouldn't experience any significant drag). Again, non-Euler (not v=f/m,s+=v) methods may improve the situation and I'll explore that later.
I've made a lot of mention of the spring and damper constants Ks and Kd, but haven't mentioned that the values must be constrained if stability will be maintained. The best thing to do is to experiment, of course, and find out what works best for the particular system. What works best is: Ks,Kd < .5If Ks or Kd were larger than this, we could have a situation where the velocity imparted will move the mass into a such a distant position in the next cycle that the spring will respond by pulling back with an even larger force, giving a velocity back that will extend the spring even further in the next update, and we get the chaos I mentioned above. It's also generally best to normalize these constants to the size of the masses. What worked decently for me for a simple sphereical arrangement was to make Kd slightly larger than Ks, while still less than .5. As it stands, my system has a uniform time step for each update regardless of how long each cycle actually took, while a more extended system would measure the milliseconds and factor this in to the position update (s += v*t instead of my simplistic s+=v). If updates come too slowly, the effect will be to increase the spring constants, perhaps over .5, and cause problems (as if less than 1 fps wasn't problem enough)
Play with the demo I made: springSphere.zip Copyright © 2002-2003 Lucas W |
<-Sample Spring
Behaviour Video
